
%% MAIN FUNCTIOM

addpath(genpath(pwd));
clear;
close;


mcmcPar.seq = 1;                                       % single chains model at present
mcmcPar.arbc = 0.5;                                  % model selection index for addition or deletion of fracture   
mcmcPar.lamda = 0.1;                                  % model selection index for updating aperture
mcmcPar.kmin = 8;                                      % minimum number of fractures
mcmcPar.kmax = 17;                                      % maximum number of fractures
mcmcPar.sig = 0.1;                                    % standard deviation for measuremet error
mcmcPar.obs = load ('observations4')    ;             % 
mcmcPar.steps = 300;                              % iteration times
mcmcPar.innerloop = 3;                             % iteration times for aperture updating
mcmcPar.thin  = 2;                                 % store fracture and fx ever thin times
mcmcPar.update = 5;                               % store output to hard disk

%% --------- parameters for fracture generation -----------------

seisdata=load ('seis5mdense.txt');               % coordiantes of seismic events with noise
seisPar.buffer = 5;                            %  noise distance in seismic dataset
seisPar.iters= 5;                             % mc updating for fracture
seisPar.s_min=[20,25,30,115];                         % prior min strikes in two major directions
seisPar.s_max=[25,30,35,125];                         % prior max strikes in two major directions 
seisPar.loss_=100.0;                            % initial loss for loss updating (any large value)
seisPar.weight=5e5;                             % weight to define loss function of seismicity 

NumChain=60;
if isempty(gcp('nocreate'))
    parpool(NumChain);
end

parfor i = 1:NumChain
    [OUTPUT{i},LOSS_temperature{i},LOSS_seis{i}] = rjMCMC(i,mcmcPar, seisPar,seisdata);
end

